source("../src/communitySim.R")
source("../src/simCancer.R")

if(F) {
# how I got the 0.4899 figure
bob = function(qq) abs(qnorm(0.95, 0, sqrt(qq)) - log(10)/2)
optim(0.5, bob)
}

totalVar = 0.4899106
totalSD = sqrt(totalVar)
propGroupVar = 0.2

doSims=F


parameters = list(
	   geno=log(1.5), env=log(1.5), "geno:env"=log(2),
		 size=150000, probGeno=0.2, 
		 probEnv=0.3, genoMiss= 0.1,
		 envMiss = 0.1, oneEnvPerCommunity =T,
     cancerMissRate=0.15, falseCancerRate=0.00045, 
		 sdGroup= sqrt(propGroupVar) * totalSD, 
		 sdPerson=sqrt(1-propGroupVar)*totalSD, 
		 Enrollmentprobs = c(20000, 40000, 50000, 40000), 
		 agerange = c(35, 69),
		 Followup = 30,
		 Ncommunity =  50,
		 EqualCommunity=F ,
		 caseCohort=T,
		 controlProp=1
)
sig = 0.001

#populationData = getPopData("../data/On129Communities.csv")  
populationData = getPopData()  


deathData = list(
  "F"=read.table("../data/MortalityNocancerFemale.txt", header=T), 
  "M"=read.table("../data/MortalityNocancerMale.txt", header=T)
)
DiffCancer= list(
	"F"=lambdaDiffCancer(myfile="../data/IncidentFemaleDiffCancer.csv",
	 	agegroup),
   M=lambdaDiffCancer(myfile="../data/IncidentMaleDiffCancer.csv",
	 	agegroup)
	 )	
  

lambda =list(F=list(x=seq(30, 90, by=5), y= exp((-0.5)*0.585)*c(116.4,169.7,272.6,401.8,555.4,746.8,1003.2,1254.8,1528.7,1741.4,
1903.1,1950.2,1922.5) / 100000)  ,
M=list(x=seq(30,90, by=5), y=exp((-0.5)*0.585)*c(61.8,81.2,123.6,247.8,466.8,877.5,1403.7,2096.6,2598.4,2887.2,
3086.4, 3156.3, 2927.5)/ 100000) )


propGroupVarseq = c(0.2, 0.35,  0.5, 0.65, 0.8,  0.95)
sdGroup= sqrt(propGroupVarseq) * totalSD 
sdPerson=sqrt(1-propGroupVarseq) * totalSD

SnumberCommunity = c(15, 30 ,50, 80)

Scolour = c("black", "black", "grey", "black")

thecolours =  apply(col2rgb(Scolour)/255, 2, paste, collapse=",")  
thelwd=c(1,1,2,2)
thelty = c(1,2,1,3)

thepch=1 

Nsim = 10


caseCohort1 = seqPowerList(list("controlProp"=c(0.8,1,2,4,8)),
	parameters=parameters, CommonCancer=c("Colon"), 
	SimulationTime=Nsim, lambda=lambda, 
	populationData=populationData, verbose=T,deathData,  caseCohort=T) 
	

parameters$geno=log(2)
	
caseCohort2 = seqPowerList(list("controlProp"=c(0.8,1,2,4,8)),
	parameters=parameters, CommonCancer=c("Colon"), 
	SimulationTime=Nsim, lambda=lambda, 
	populationData=populationData, verbose=T,deathData,  caseCohort=T) 
	
parameters$geno=log(4)
	
caseCohort3 = seqPowerList(list("controlProp"=c(0.8,1,2,4,8)),
	parameters=parameters, CommonCancer=c("Colon"), 
	SimulationTime=Nsim, lambda=lambda, 
	populationData=populationData, verbose=T,deathData,  caseCohort=T) 
